Energy transfer in nonlinear network models of proteins 
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We investigate how nonlinearity and topological disorder affect the energy relaxation of local 
kicks in coarse-grained network models of proteins. We find that nonlinearity promotes long-range, 
coherent transfer of substantial energy to specific, functional sites, while depressing transfer to 
generic locations. Remarkably, transfer can be mediated by the self-localization of discrete breathers 
at distant locations from the kick, acting as efficient energy-accumulating centers. 
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It is now well established that the functional dynam- 
ics of proteins is deeply rooted in the peculiar topo- 
logical arrangement of their native folds, as revealed 
by many experimental and computational studies [TJ|2]. 
In particular, the success of coarse-grained elastic net- 
work models (ENMs) in describing atomic fluctuations 
at room temperature have helped elucidate, at the har- 
monic level, the subtle interplay between structure and 
dynamics on one side and biological function on the 

other |3H00HEig|TDj. 

However, protein dynamics is strongly anharmonic p~2j 
[13], a property which has to be taken into account in 
order to rationalize crucial biological processes such as 
energy storage and transfer upon ligand binding, chem- 
ical reaction, etc [HI [15] • Yet, even though many the- 
oretical studies suggest that nonlinear excitations may 
play an active role in protein functioning [16] [17], the 
rich phenomenology residing in the interplay between 
protein topology and nonlinearity still remains widely 
unexplored. Along these lines, we have recently intro- 
duced the Nonlinear Network Model (NNM), showing 
how known nonlinear effects can be modulated by the 
underlying non-regular topology of protein systems. For 
instance, within a large collection of enzyme structures, 
the formation of localized, robust nonlinear modes ap- 
pears strongly favored at few specific sites, that often lie 
in close proximity of known catalytic sites [T8l [T9], 

In this paper we examine the effects of the nonlin- 
earity /topology interplay on energy transfer phenomena 
across protein structures. Within the NNM framework a 
protein is represented by N fictive particles (amino acids) 
of identical mass M = 110 a.m.u., at equilibrium at the 
corresponding C a site as specified in the experimentally 
determined structures (X-ray or NMR). By imposing a 
fixed cutoff R c on the latter set of coordinates, a protein 
is mapped onto a network of nonlinear oscillators, whose 
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where is the distance between residues i and j, R^ 
their distance in the equilibrium structure and Cij = 
{1 if Rij < i? c ,0 otherwise} is the connectivity matrix. 
As in previous studies [19], we take R c =10 A, k± = 5 
kcal/mol/A 4 and fix so that the low- frequency part of 
the linear spectrum match actual protein frequencies, as 
calculated through realistic force fields [20] [2TJ [22] . This 
gives &2 = 5 kcal/mol/A 2 . The case k± = corresponds 
to the Anisotropic Network Model (ANM) [3] [4] [5]. 

Our aim is to investigate how energy initially imparted 
at a specified site i redistributes across a given structure. 
To do this, we perform microcanonical simulations with 
all residues initially at rest at their equilibrium position 
but for a kinetic energy kick at site i of magnitude E . 
Sites in a 3D protein network are not equivalent, featur- 
ing e.g. varying connectivity, clustering coefficient and 
bond directions. Thus, in order to allow for a compari- 
son of energy relaxation from all sites in a given structure, 
the initial kick direction ought to be specified by a unique 
protocol. We chose to calculate the directions of the 
initial velocities Vi(0) through the Sequential Maximum 
Strain (SMS) algorithm [18], which provides an unbiased 
measure of the maximum-strain direction at site i for a 
fixed displacement (here 1 A), e SMS . During the simu- 
lation, we record at regular intervals k = 1, 2, . . . , N s 
the site that carries the highest energy, n/e, and the value 
of the latter, e Uk = Mv^ k /2 + u nk . Corresponding to a 
fixed simulation time (about 500 ps), we define a transfer 
probability from site i (the kicked one) to site j and the 
fraction of energy transferred as 
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Figure 1: Average transfer probability in Riboflavin Synthase 
(PDB id. 1KZL, N = 202) on logarithmic (upper panel) and 
linear (lower panel) scale. Thick solid line: ANM. Dashed 
line: NNM. The staircase plot in the lower panel reproduces 
the stiffness pattern (arb. units). Filled circles flag catalytic 
sites. Eq — 75 kcal/mol. 



The first striking result comes from the calculation of 
average transfer probabilities. These gauge the mean 
transfer to a given site from kicks at all other sites, 
(Vj) = Vi^j / (N — 1), obtained from N indepen- 

dent simulations. A typical probability transfer plot is 
shown in Fig. [I] The first notable feature is that the ef- 
fect of nonlinearity is to substantially increase the prob- 
ability of energy funneling to a few selected sites, while 
depressing transfer to all other locations with respect to 
the harmonic (ANM) case. Remarkably, the preferred 
target sites lie in close proximity to the known catalytic 
sites, within the stiffest regions [26]. Thus, topology and 
nonlinearity team in this case together to sharpen energy 
funneling to specific functional regions. 

The case shown in Fig.[l]is not a singular one. In Fig. [2] 
we show the stiffness patterns for four other enzymes 
along with the sites ranking first to tenth as to the energy 
delivered on average to spherical shells with 6 A radius 
around each site. For residue j, this amounts to further 
averaging the mean energy deposited at sites within the 
j-th ball B(j), i.e. «//» = (Ei&f^j/iN - l)) m . 
It is manifest that the sites around which most of the 
energy is deposited invariably spotlight the stiffest re- 
gions, at the same time identifying functionally relevant 
locations (see catalytic sites). Moreover, the same lo- 
cations clearly attract substantial fractions of the initial 
excitation energy, as revealed by surveying the maximum 
transferred energies to each ball B(j), that is / m ax(j) = 
(max^j fi-^j)B(j) (empty circles). Many events featur- 
ing transfers of energy fractions in the range 20 to 25 % 



were indeed observed. 

We can learn more on the mechanisms underlying the 
energy transfer process by examining in detail the out- 
come of a single kick. Fig. [3] pictures a long-range trans- 
fer event occurring when kicking at site LEU 42 in the 
enzyme Subtilisin. The middle-lower panel (a) shows a 
plot of the most energetic site as a function of time, 
clearly illustrating the transfer to site VAL 177, some 
23 A away, occurring at t* « 275 ps. The transfer pro- 
cess also involves site ALA 85, as a passage site. Re- 
markably, a plot of the energy e max (t) of the most ener- 
getic site at time t clearly shows that such passage co- 
incides with a redistribution of energy across the struc- 
ture (see middle panel (c)). Subsequently, energy is gar- 
nered from the neighborhood and stabilized in a local- 
ized mode centered at VAL 177, finally carrying about 
20 % of the total energy. This marks the true transfer 
event. Such energy-harvesting, self- localized vibrations 
are generic in discrete non-linear systems and are well- 
known as Discrete Breathers (DB) [23]. These are ro- 
bust, time-periodic exponentially localized modes, whose 
vibrational frequency lies outside the linear spectrum of 
the system. In the context of the NNM, we have shown 
how accurate approximations of such periodic orbits can 
be calculated analytically, reproducing the marked affin- 
ity of DB self-localization in topologically disordered me- 
dia for the stiffest spots [181 [19]. Here we have shown 
that DBs may also be excited as a consequence of lo- 
calized impulses at considerable distances from the exci- 
tation, playing the role of energy-accumulating transfer 
vectors. In order to substantiate the above interpreta- 
tion, we have performed Principal Component Analysis 
(PCA) on an extended portion of the post-transfer dy- 
namics. The power spectrum of the system trajectory 
projected on the first principal mode (PM1) is shown 
in the upper left panel of Fig. [3J clearly revealing the 
existence of a nonlinear, time-periodic excitation, remi- 
niscent of chaotic DBs that self-localize through modu- 
lational instability in nonlinear lattices [24] . 

More insight as to why energy is transferred from LEU 
42 a long distance away can be obtained by turning to the 
space of Normal Modes (NM) ^ (k = 1, 2, . . . , 3N - 6). 
Middle panel (b) of Fig. p] reports the mode carrying the 
highest energy = (Q| +ci;^Q|)/2 at time t, where 

Qk = a %ja(,j a i s the NM- transform of the system 
coordinates Xj a (j = 1, 2, . . . , N, a = x, y : z) and uok are 
the NM frequencies. Before the transfer event, energy 
is bounced among four high-frequency modes, NMs 1,3,8 
and 9. This can be understood by constructing the NM 
overlap network starting from the NMs with the largest 
projections on the initial excitation unit vector e SMS are 
greatest (four modes in the first column of the bottom 
graph in Fig.[3] NMs 138,127,67 and 37, making up about 
60 % of e SMS ). For a given NM p, two links are drawn 
to the two-highest ranking NMs in the ordered list of 
absolute overlap coefficients t pq = Y^j, a IfjUlfjU- B ^ 
doing this for the four NMs involved in the SMS vector, 
a closed network emerges identifying the NM3 ^ NM8 ^ 
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Figure 2: Stiffness plots (staircases) and first ten sites in the ranked list of average energies delivered to 6 A-balls around each 
site (dashed impulses) for four enzymes: HIV-I protease (PDB id. 1A30, N = 201), Astacin (PDB id. 1AST, N = 200), 
Imidazole glycerol phosphate synthase subunit hisF (PDB id. 2A0N, N = 200), SARS coronavirus main proteinase (PDB id. 
2BX4, N = 299). Filled circles flag catalytic sites. Right axes report the ten largest fractions of energy transferred (in units of 
.Eo), averaged over 6 A balls, versus target ball centers (empty circles). Eq = 75 kcal/mol. 



NM9 loop. Thus, in the presence of nonlinearity energy 
is immediately directed to a reduced group of NMs via 
resonant overlap mechanisms. This finding agrees with 
results of atomistic simulations highlighting the impor- 
tance of spatial overlap for NM-NM energy transfer [25] . 

High-frequency NMs are strongly localized in space. In 
particular, ALA 85 is the NM site (the site with largest 
displacement) in NM3 and the second NM site in NM8, 
which explains the role of ALA 85 in the energy cir- 
culation process. Before transfer, however, energy also 
bounces back and forth from NM1, the highest-frequency 
mode, reflecting the nonlinear frequency shift on NM3 to- 
ward greater frequencies (see again panel b). At t = t* 
energy starts departing the region around LEU 42 and a 
fluctuation pumping up NM3 occurs (panel d), shifting 
its frequency upwards by virtue of nonlinearity. The en- 
ergy at stake is sufficient to trigger nonlinear localization 
and a DB finally installs at VAL 177, the NM site of NM1, 
gathering vibrational energy from the background. Cor- 
respondingly, the energy on NM1 increases (see panel d). 
To substantiate the above analysis, we have calculated 
analytically the DB mode pattern centered at site VAL 
177 with the technique described in Ref. [18| Then we 
have built the network connecting the first two principal 
modes, the first three NMs and the DB, where the links 
are weighted by the normalized scalar products (upper 



graph in Fig.[3|. As it shows, the PMs essentially reflect 
the underlying competition between NM1 and NM3. In 
particular, the first principal mode confirms the excita- 
tion of a DB emerging as a nonlinear continuation of the 
edge normal mode, as predicted theoretically in Ref. [THl 
In agreement with this picture, kicks at ALA 42 of weaker 
energy resulted in a DB installing at MET 199, the NM 
site of NM2. That is, less energy causes a smaller fre- 
quency shift and the DB branch originating from the 
continuation of NM2 is excited instead. Reducing Eq 
further, the transfer is observed to halt at ALA 85, as 
explained by the NM overlap network. 

In this paper we have shown how nonlinearity in a 
topologically non-regular system boosts energy transfer 
to few specific locations. In enzyme structures, these 
coincide invariably with the stiffest regions, also host- 
ing the functionally relevant sites. Nonlinearity sharp- 
ens the transfer selectivity, by reducing at the same time 
the transfer probability to generic locations. The energy 
transferred by virtue of nonlinearity may be a conspic- 
uous portion of the initial excitation, in which cases lo- 
calized vibrations akin to Discrete Breathers self-localize 
as energy-collecting centers, often realizing amazingly ef- 
ficient energy transfer channels across considerable dis- 
tances. 
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Figure 3: (Color online) Kick at site LEU 42 in Subtilisin (PDB id. 1AV7, N = 274). Middle panels: most energetic site vs 
time (a), most energetic normal mode vs time - NMs being ranked in order of decreasing eigenfrequency - (b), highest site 
energy vs time (c), energies of two NMs vs time (d). Upper left panel: power spectrum of the system trajectory for t > £* = 275 
ps projected on the first principal mode; uj m — 101.2 cm -1 is the band-edge linear frequency. Lower right graph: NM overlap 
network. Nodes are NMs, red and blue links connect to nodes ranking first and second, respectively, in overlap (see text). Link 
weights are the overlap coefficients t pq . Upper right graph. Network relating principal modes (PM), NMs and the analytical 
Discrete Breather pattern (DB). Link weights are the absolute cosines (normalized scalar products). In both graphs the link 
width is proportional to its weight. Eq = 100 kcal/mol. 
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